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Abstract 

We investigate the detailed dynamics of multidimensional Hamiltonian systems by 
studying the evolution of volume elements formed by unit deviation vectors about 
their orbits. The behavior of these volumes is strongly influenced by the regu- 
lar or chaotic nature of the motion, the number of deviation vectors, their linear 
(in)dependence and the spectrum of Lyapunov exponents. The different time evolu- 
tion of these volumes can be used to identify rapidly and efficiently the nature of the 
dynamics, leading to the introduction of quantities that clearly distinguish between 
chaotic behavior and quasiperiodic motion on A^-dimensional tori. More specifically 
we introduce the Generalized Alignment Index of order k (GALI^) as the volume of 
a generalized parallelepiped, whose edges are k initially linearly independent unit 
deviation vectors from the studied orbit whose magnitude is normalized to unity at 
every time step. We show analytically and verify numerically on particular examples 
of N degree of freedom Hamiltonian systems that, for chaotic orbits, GALI^ tends 
exponentially to zero with exponents that involve the values of several Lyapunov 
exponents. In the case of regular orbits, GALI^ fiuctuates around non-zero values 
for 2 < k < N and goes to zero for N < k < 2N following power laws that depend 
on the dimension of the torus and the number m of deviation vectors initially tan- 
gent to the torus: oc t-'^{k-N)+m if q < m < A; - A^, and tx ^-C^-^) iim> k- N. 
The GALIfc is a generalization of the Smaller Alignment Index (SALI) (GALI2 oc 
SALI). However, GALI^ provides significantly more detailed information on the lo- 
cal dynamics, allows for a faster and clearer distinction between order and chaos 
than SALI and works even in cases where the SALI method is inconclusive. 
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1 Introduction 



Determining the chaotic or regular nature of orbits in conservative dynamical 
systems is a fundamental issue of nonlinear science. The difficulty with con- 
servative systems, of course, is that regular and chaotic orbits are distributed 
throughout phase space in very complicated ways, in contrast with dissipa- 
tive systems, where all orbits eventually fall on regular or chaotic attractors. 
Over the years, several methods distinguishing regular from chaotic motion in 
conservative systems have been proposed and applied, with varying degrees 
of success. These methods can be divided in two major categories: Some are 
based on the study of the evolution of small deviation vectors from a given 
orbit, while others rely on the analysis of the particular orbit itself. 

The most commonly employed method for distinguishing between order and 
chaos, which belongs to the category related to the study of deviation vectors, 
is the evaluation of the maximal Lyapunov Characteristic Exponent (LCE) ai; 
if 0"! > the orbit is chaotic. The theory of Lyapunov exponents was applied 
to characterize chaotic orbits by Oseledec [T], while the connection between 
Lyapunov exponents and exponential divergence of nearby orbits was given in 
[2][3] . Benettin et al. [1] studied the problem of the computation of all LCEs 
theoretically and proposed in [5] an algorithm for their numerical computation. 
In particular, cxi is computed as the limit for t — >^ cx) of the quantity 



where w{0), w{t) are deviation vectors from a given orbit, at times t = and 
t > respectively. It has been shown that the above limit is finite, independent 
of the choice of the metric for the phase space and converges to cxi for almost all 
initial vectors w{0) [Tf^fS] . Similarly, all other LCEs, cr2, us etc. are computed 
as the limits for t — > cx) of some appropriate quantities, L2(t), L^{t) etc. (see 
[5] for more details). We note that throughout the present paper, whenever we 
need to compute the values of the maximal LCE or of several LCEs we apply 
respectively the algorithms proposed by Benettin et al. [211^. Since 1980, new 
methods have been introduced for the effective computation of LCEs (e. g. 
[6], see also [7] and references therein). The true power of these techniques 
is revealed in the study of multi-dimensional systems, when only a small 
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number of LCE are of interest. In such cases, these methods are significantly 
more efficient than the method of [5j, which computes the whole spectrum of 
LCEs. On the other hand, they are less or equally efficient when compared 
with the method of [2j for the computation of the maximal LCE, whose value 
is sufficient for the determination of the regular or chaotic nature of an orbit. 

Among other chaoticity detectors, belonging to the same category with the 
evaluation of the maximal LCE, are the fast Lyapunov indicator (FLI) and its 
variants Pl9lll0|lllp2| . the mean exponential growth of nearby orbits (MEGNO) 
[131111], the smaller alignment index (SALI) [TSlll6fl7j . the relative Lyapunov 
indicator (RLI) [TS], as well as methods based on the study of power spec- 
tra of deviation vectors |19], as well as spectra of quantities related to these 
vectors [20|21|22j . In the category of methods based on the analysis of a time 
series constructed by the coordinates of the orbit under study, one may list 
the frequency map analysis of Laskar |23|24|25|26II27|28] . the method of the 
low frequency power (LFP) [2^1150] . the '0-1' test [HI], as well as some other 
more recently introduced techniques [321133] . 

In the present paper, we generalize and improve considerably the SALI method 
mentioned above by introducing the Generalized ALignment Index (GALI). 
This index retains the advantages of the SALI - i.e. its simplicity and effi- 
ciency in distinguishing between regular and chaotic motion - but, in addi- 
tion, is faster than the SALI, displays power law decays that depend on torus 
dimensionality and can also be applied successfully to cases where the SALI 
is inconclusive, like in the case of chaotic orbits whose two largest Lyapunov 
exponents are equal or almost equal. 

For the computation of the GALI we use information from the evolution of 
more than two deviation vectors from the reference orbit, while SALI's com- 
putation requires the evolution of only two such vectors. In particular, GALI^ 
is proportional to 'volume' elements formed by k initially linearly independent 
unit deviation vectors whose magnitude is normalized to unity at every time 
step. If the orbit is chaotic, GALI^ goes to zero exponentially fast by the law 



If, on the other hand, the orbit lies in an A^-dimensional torus, GALI^ displays 
the following behaviors: Either 



or, if < A; < 2A^, it decays with different power laws, depending on the 
number m of deviation vectors which initially lie in the tangent space of the 



GALIfc(t) oc e 



[(o-l-o-2) + (o"l-o-3)H \-io-i-a-k)]t 



(2) 



GALIfc(t) ^ constant for 2 < k < N, 



(3) 
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torus, i. e. : 



GALIfc(t) oc <^ * ' ' (4) 
^ if N < k < 2N and m > k - N 



So, the GALI allows us to study more efficiently the geometrical properties of 
the dynamics in the neighborhood of an orbit, especially in higher dimensions, 
where it allows for a much faster determination of its chaotic nature, overcom- 
ing the limitations of the SALI method. In the case of regular motion, GALI^ 
is either a constant, or decays by power laws that depend on the dimensional- 
ity of the subspace in which the orbit lies, which can prove useful e.g., if our 
orbits are in a 'sticky' region, or if our system happens to possess fewer or 
more than N independent integrals of the motion (i.e. is partially integrable 
or super- integrable respectively). 

This paper is organized as follows: In section [2|, we recall the definition of 
the SALI describing also its behavior for regular and chaotic orbits of Hamil- 
tonian flows and symplectic maps. In section [3], we introduce the GALI^ for 
k deviation vectors, explaining in detail its numerical computation, while in 
section H] we study theoretically the behavior of the new index for chaotic 
and regular orbits. Section [5] presents applications of the GALIfe approach to 
various Hamiltonian systems of different numbers of degrees of freedom, con- 
centrating on its particular advantages. Finally, in section [6l we summarize 
the results and present our conclusions, while the appendices are devoted re- 
spectively to the definition of the wedge product and the explanation of the 
explicit connection between GALI2 and SALI. 



2 The SALI 



The SALI method was introduced in JT5\ and has been applied successfully to 
detect regular and chaotic motion in Hamiltonian flows as well as symplectic 
maps [M|16ll35|36|17|37|38|39ll40f41f42ii43.44j . It is an index that tends ex- 
ponentially to zero in the case of chaotic orbits, while it fluctuates around 
non-zero values for regular trajectories of Hamiltonian systems and 2N- 
dimensional symplectic maps with > 1. In the case of 2-dimensional (2D) 
maps, the SALI tends to zero both for regular and chaotic orbits but with 
very different time rates, which allows us again to distinguish between the two 
cases [15]: In particular the SALI tends to zero following an exponential law 
for chaotic orbits and decays to zero following a power law for regular orbits. 

The basic idea behind the success of the SALI method (which essentially 
distinguishes it from the computation of LCEs) is the introduction of one 
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additional deviation vector with respect to a reference orbit. Indeed, by con- 
sidering the relation between two deviation vectors (instead of one deviation 
vector and the reference orbit), one is able to circumvent the difficulty of 
the slow convergence of Lyapunov exponents to non-zero (or zero) values as 
t ^ oo. 

In order to compute the SALI, therefore, one follows simultaneously the time 
evolution of a reference orbit along with two deviation vectors with initial 
conditions wi(0), W2(0). Since we are only interested in the directions of these 
two vectors we normalize them, at every time step, keeping their norm equal 
to 1, setting 

Mt) = T^,, ^ = l,2 (5) 



where || ■ || denotes the Euclidean norm and the hat (^) over a vector denotes 
that it is of unit magnitude. The SALI is then defined as: 

SALI(t) = min{||^i(t) + W2{t)\\ , \\wi{t) - Mt)\\} , (6) 



whence it is evident that SALI(t) G [0, \/2]. SALI = indicates that the two 
deviation vectors have become aligned in the same direction (and are equal or 
opposite to each other); in other words, they are linearly dependent. 

Let us observe, at this point, that seeking the minimum of the two positive 
quantities in ([6]) (which are bounded above by 2) is essentially equivalent to 
evaluating the product 

Pit) = \\wi{t) + W2{t) II ■ \\Wi{t) - W2{t) II , (7) 



at every value of t. Indeed, if the minimum of these two quantities is zero (as 
in the case of a chaotic reference orbit, see below), so will be the value of P{t). 
On the other hand, if it is not zero, P{t) will be proportional to the constant 
about which this minimum oscillates (as in the case of regular motion, see 
below). This suggests that, instead of computing the SALI(t) from ([6]), one 
might as well evaluate the 'exterior' or 'wedge' product of the two deviation 
vectors Wi A W2 for which it holds 

II , II ll^^i - ^all ■ + W2\\ 

\\Wi /\W2\\ = , (8) 



and which represents the 'area' of the parallelogram formed by the two devia- 
tion vectors. For the definition of the wedge product see Appendix |A] and for 
a proof of (IE]) see Appendix [Bl Indeed, the 'wedge' product can provide much 
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more useful information, as it can be generalized to represent the 'volume' of 
a parallelepiped formed by the vectors Wi,W2,---,Wk, 2 < k < 2N, regarded 
as deviations from an orbit of an A^-degree of freedom Hamiltonian system, 
or a 2iV-dimensional symplectic map. 

It is the main purpose of this paper to study precisely such a generalization 
and reveal considerably more qualitative and quantitative information about 
the local and global dynamics of these systems. Before we proceed to describe 
this generalization, however, let us first summarize what we know about the 
properties of the SALI for the case of two deviation vectors Wi, W2: 

(1) In the case of chaotic orbits, the deviation vectors wi, W2 eventually be- 
come aligned in the direction of the maximal Lyapunov exponent, and 
SALI(t) falls exponentially to zero. An analytical study of SALI's behav- 
ior for chaotic orbits was carried out in [T7| where it was shown that 

SALI(t) oc e-("^-"2)* (9) 

(Ji, (72 being the two largest LCEs. 

(2) In the case of regular motion, on the other hand, the orbit lies on a torus 
and the vectors wi, W2 eventually fall on its tangent space, following a 

time evolution, having in general different directions. In this case, the 
SALI oscillates about values that are different from zero (for more details 
see [IS]). This behavior is due to the fact that for regular orbits the norm 
of a deviation vector increases linearly in time along the flow. Thus, our 
normalization procedure brings about a decrease of the magnitude of the 
coordinates perpendicular to the torus at a rate proportional to and 
so wi, W2 eventually fall on the tangent space of the torus. 

Note that in the case of 2D maps the torus is actually an invariant curve and 
its tangent space is 1-dimensional. So, in this case, the two unit deviation vec- 
tors eventually become linearly dependent and SALI becomes zero following 
a power law. This is, of course, different than the exponential decay of SALI 
for chaotic orbits and thus SALI can distinguish easily between the two cases 
even in 2D maps [15j. Thus, although the behavior of SALI in 2D maps is 
clearly understood, the fact remains that SALI does not always have the same 
behavior for regular orbits, as it may oscillate about a constant or decay to 
zero by a power law, depending on the dimensionality of the tangent space 
of the reference orbit. It might, therefore, be interesting to ask whether this 
index can be generalized, so that different power laws may be found to char- 
acterize regular motion in higher dimensions. It is one of the principal aims of 
this paper to show that such a generalization is possible. 

Let us make one final remark concerning the behavior of SALI for chaotic 
orbits: Looking at equation IQ, one might wonder what would happen in 
the case of a chaotic orbit whose two largest Lyapunov exponents cti and 
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(72 are equal or almost equal. Although this may not be common in generic 
Hamiltonian systems, such cases can be found in the literature. In one such 
example [39], very close to a particular unstable periodic orbit of a 15 degree of 
freedom Hamiltonian system, the two largest Lyapunov exponents are nearly 
equal cti — o"2 ~ 0.0002. Even though, in that example, SALI still tends to zero 
at the rate indicated by (Q, it is evident that the chaotic nature of an orbit 
cannot be revealed very fast by the SALI method. It is, therefore, clear that 
a more detailed analysis of the local dynamics is needed to further explore 
the properties of specific orbits, remedy the drawbacks and improve upon 
the advantages of the SALI. For example, if we could define an index that 
depends on several Lyapunov exponents, this might accelerate considerably 
the identification of chaotic motion. 



3 Definition of the GALI 



Let us consider an autonomous Hamiltonian system of degrees of freedom 
having a Hamiltonian function 

H{qi,q2, . . .,qN,Pi,P2, . . . ,Pn) = h = constant (10) 



where Qi and pi, i = 1,2, N are the generalized coordinates and conjugate 
momenta respectively. An orbit of this system is defined by a vector x{t) = 

{qi{t),q2{t),---,qNit),Pi{t),p2{t),...,PN{t)), with Xi = Qi, Xi+N = Pi, i = 

1,2, . . . , N. The time evolution of this orbit is governed by Hamilton equations 
of motion 

dx ^ (dH dH\ 



while the time evolution of an initial deviation vector w(0) = {dxi{0), . . . , c?X2Ar(0)) 
from the x{t) solution of ffTTj) obeys the variational equations 

^ = M(x(t))iU, (12) 



where M = dV /dx is the Jacobian matrix of V. 

The SALI is a quantity suitable for checking whether or not two normalized 
deviation vectors Wi, W2 (having norm 1), eventually become linearly depen- 
dent, by falling in the same direction. The linear dependence of the two vectors 
is equivalent to the vanishing of the 'area' of the parallelogram having as edges 
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the two vectors. Generalizing this idea we now follow the evolution of k de- 
viation vectors Wi, W2, ■ ■ ■, w^, with 2 < k < 2N, and determine whether 
these eventually become linearly dependent, by checking if the 'volume' of the 
parallelepiped having these vectors as edges goes to zero. This volume will be 
computed as the norm of the wedge product of these vectors (see Appendix 
|X]for a definition of the wedge product). 

All normalized deviation vectors Wi, i = 1,2, . . . , k, belong to the 2iV-dimensional 
tangent space of the Hamiltonian flow. Using as a basis of this space the usual 
set of orthonormal vectors 

ei = (l,0,0,...,0),e2 = (0,l,0,...,0),...,e2Ar = (0,0,0,...,!) (13) 
any deviation vector Wi can be written as 

2N 

'^i = J2'^ij^j ' ^ = l,2,...,/c (14) 
i=i 

where Wij are real numbers satisfying 

2N 

Thus, equation (1A.12P of Appendix \M gives 



Wi 




Wii W12 ■ 


■ Wi2N 




ei 




W2 




W21 W22 ■ 


■ W22N 






(16) 


Wk 




Wkl Wk2 ■ 


■ Wk2N 




e2N 





Using then equation (1A.13P the wedge product of these k deviation vectors 
takes the form 



Wi A t&2 A ■ ■ ■ A t&fc = 

l<'ti<«2<--<ifc<2Af 



Wii^ Wii2 ■ ■ ■ Wii^ 
W2h ^2i2 ■ ■ ■ W2i^ 

Wkn Wki2 ■ ■ ■ Wki^ 



Ci, A A • ■ ■ A Ci^ , (17) 
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where the sum is performed over all possible combinations of k indices out of 
2N. 



If at least two of the normalized deviation vectors Wi,i = 1,2, ... ,k are linearly 
dependent, all the k x k determinants appearing in equation (IT7|) will become 
zero making the 'volume' vanish. Equivalently the quantity 



\wi A W2 A ■ ■ ■ A Wk\ 



E 

l<ii<i2<---<ife<2Af 



Wii^ Wii^ ■ ■ ■ Wu^ 



Wk 



2^ 



1/2 



(18) 



which we shall call the 'norm' of the wedge product, will also become zero. 
Thus, we define this important quantity as the Generalized Alignment Index 
(GALI) of order k 

GALIfc(t) = AwsW A ■■■ AWfc(t)|| . (19) 



In order to compute GALI^, therefore, we need to follow the evolution of an 
orbit with initial conditions ^(0), using equation ffTTj) . as well as the evolution 
of k initially linearly independent unit deviation vectors Wi, i = 1,2, ... ,k 
using the variational equations f[T^ . At every time step, we normalize these 
deviation vectors to unity and compute GALI^ as the norm of their wedge 
product using equation (fT8|) . 

Consequently, if GALlkit) tends to zero, this would imply that the volume of 
the parallelepiped having the vectors Wi as edges also shrinks to zero, as at least 
one of the deviation vectors becomes linearly dependent on the remaining ones. 
On the other hand, if GALIfc(t) remains far from zero, as t grows arbitrarily, 
this would indicate the linear independence of the deviation vectors and the 
existence of a corresponding parallelepiped, whose volume is different from 
zero for all time. 



4 Theoretical results 



4-1 Exponential decay of GALI for chaotic orbits 



In order to investigate the dynamics in the vicinity of a chaotic orbit of 
the Hamiltonian system (llOp with N degrees of freedom, let us first recall 
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some known properties of the Lyapunov characteristic exponents, following e. 
g. [ISPG] . It has been shown that the mean exponential rate of divergence 
a (x(0), iS) from a reference orbit with initial condition x(0) given by 



a(f(0),^)= limllnMziL, (20) 



exists and is finite. Furthermore there is a 2A^-dimensional basis {u\, U2, ■ ■ ■ , U2n} 
of the tangent space of the Hamiltonian flow so that a {x{0), tv) takes one of 
the 2N (possibly nondistinct) values 

ai{x{0)) = a{x{0),Ui) , t = l,2,...,2N (21) 



which are the Lyapunov characteristic exponents, ordered in size as follows: 

ai > (T2 > . . . > . (22) 



These properties can be easily understood if the reference orbit is an unstable 
periodic solution of period T. In this case, the matrix M of the variational 
equations (fT2!) is a continuous T-periodic 2N x 2N matrix. The solution of 
equations (fT2l) can be written as 

w{t) = $(t) ■ w{0) , (23) 



where $(t) is the so-called fundamental matrix (see e. g. [S]), such that 
$(0) = I, the 2N X 2N identity matrix. The behavior of the deviation vec- 
tor w{t) and consequently the stability of the periodic orbit is determined 
by the eigenvalues A, of the so-called monodromy matrix $(T), ordered as 
|Ai| > IA2I > •■■ > |A2Ar|. Let Ui, i = 1,2, . . . ,2N denote the corresponding 
eigenvectors. Then for w{0) = Ui we have 

w{nT) = X^Ui , i = l,2,...,2N (24) 

and from (!20l) we get 

a(f(0),M,) = lim4pln|Ar| = ^ , i = 1,2, . . . ,2N. (25) 

Furthermore, if we write 

2N 

wiO) =Y.CiUi, (26) 
1=1 
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it follows from that the first nonvanishing coefficient q dominates the 
subsequent evolution of w{nT). Thus, if Ci 7^ we get from (!20|) a (x(0), w) = 
(Ti, if Ci = and C2 7^ we get a (x(0), w) = a2 and so on. So, the evolution 
of the initial deviation vector w{0) is well approximated by 



2N 



W 



[uT) 



(27) 



For a nonperiodic orbit we cannot define such eigenvalues and eigenvectors 
as above. Nevertheless, Oseledec has proven the existence of basis vec- 
tors {-ui, ^2, . . . , U2n} and Lyapunov exponents for nonperiodic orbits. This is 
perhaps not surprising, since periodic orbits are dense in the phase space of 
Hamiltonian systems and thus a periodic orbit of arbitrary large period can 
always be found arbitrary close to any nonperiodic orbit. So, the time evolu- 
tion of a deviation vector may be approximated by a variant of equation fl27p . 
i. e. 



w{t) 



2N 
i=l 



where Cj, di are real numbers depending on the specific phase space loca- 
tion through which the reference orbit passes. Thus, the quantities di, i = 
1,2,..., 2N may be thought of as 'local Lyapunov exponents' having as limits 
for t — s> 00 the 'global' LCEs at, i = 1,2,..., 2N. We notice that even if in 
some special cases where the vectors Ui, i = 1,2, ... , 2N are known a priori, 
so that one could set w{0) = Ui, the computational errors in the numerical 
evolution of the deviation vector would lead to the actual computation of cti 
from equation (JT]) [5j. 

It is well known that Hamiltonian systems are generically non-integrable and 
possess Lyapunov exponents in chaotic domains which are real and grouped in 
pairs of opposite sign with two of them being equal to zero. We, therefore, have 
o"i = -(72N-i+i for i = 1, 2, . . . , and cti > 0-2 > ■ ■ ■ > cTAf-i >(^n = <^n+i = 

> crN+2 > ■ • • > cr2N- Assuming that, after a certain time interval, the di, 

1 = 1,2, . . . , 2N do not fiuctuate significantly about their limiting values, we 
write di ~ (jj and express the evolution of the deviation vectors Wi in the form 

2N 

w^{t) = Y^c]e''^'u, (29) 



(see discussion in section |5TT] and figured]). Thus, if ai > a2, a leading order 
estimate of the deviation vector's Euclidean norm (for t large enough), is given 
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by 



(30) 



Consequently, the matrix C in (]A.12p of coefficients of k normalized deviation 
vectors Wiit) = Wi(t)/\\wi(t)\\, i = 1,2, . . . , k with 2 < k < 2N, using as basis 
of the vector space the set {iti, U2, ■ ■ ■ , U2n} becomes 



C{t) 



Sl 



3_^ — (o-i— CT2)t (0-1—0-3)4 £2JV^-(oi— 02jv)i 



\c\\ 



2 ^ i 

^2 J^2_g— (o"i— 02)* _£^g-(oi— 03)4 . . . £2W.g-(oi— 02jv)i 



14 
^2 



,2 



i_p-(o"l-o"2)t _£3 ^-(01-03)4 . . . !hNp-{a-i-a2N)t 



\ct\ 



(31) 



with Si = sign(c*]^) and i = 1,2, . . . , k, j = 1,2, . . . , 2N and so we have 



Wi W2 . . . Wk 



C 



U1U2 . . . U2N 



(32) 



with C") denoting the transpose of a matrix. The wedge product of the k 
normalized deviation vectors is then computed as in equation f |T71) by: 



wi{t) Aw2it) A ■ • • A Wkit) 



E 

l<n<i2<---<ifc<27V 



Mj, A A • • ■ A Mifc.(33) 



Note that the quantity 



Sk={ 



E 

l<ii<j2<--<ifc<2iV 



C2n C2i2 



2^ 



1/2 



(34) 



is not identical to the norm (ITS!) of the /c-vector i£'i()f:) A W2(i':) A ■ ■ ■ A ^^(/f:) as 
the wedge product in equation (15^ is not expressed with respect to the basis 
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f ll3p . Thus one should consider the transformation 



U1U2 . . . U2N 



ei 62 



e2N 



(35) 



between the two bases, with Tc denoting the transformation matrix. Of course, 
when considering the wedge product of 2N deviation vectors one can easily 
show that 



\wi A A ■ • • A W2n\\ = S'_ 



2N 



IdetT, 



(36) 



If, on the other hand, we consider the wedge product of fewer than 2N devia- 
tion vectors, the norm (ITS!) and the quantity 5*^ are not related through a 
simple expression like f l36l) . We shall proceed, however, to obtain results using 
f l34l) instead of (fTSl) . as we do not expect that such a change of basis will affect 
significantly the dynamics and alter our conclusions for the following reasons: 
First, we note that both quantities are zero when at least two of the k devia- 
tion vectors are linearly dependent, due to the fact that all the determinants 
appearing in equations f[T5]) and fl5^ vanish. In addition, the transformation 
matrix is not singular as the sets {iii} and {cj}, 2 = 1,2,..., 2N continue 
to be valid bases of the vector space. Thus, both quantities are expected to 
behave in a similar way in the case of chaotic orbits, where the deviation vec- 
tors tend to become linearly dependent. Thus, by studying analytically the 
time evolution of Sk through (IMI) . we expect to derive accurate approxima- 
tions of the behavior of the GALI^ (1191) for chaotic orbits. The validity of this 
approximation is numerically tested and verified in section [51 



Let us now see how this approximation is derived: The determinants appearing 
in the definition of Sk (see equation (IMI) ) can be divided in two categories 
depending on whether or not they contain the first column of matrix C. Using 
standard properties of determinants, we see that those that do contain the 
first column yield 



\cf\ 



Jfc-ig-{'^i-'^jfe_i)« 
^1 ^ 



\ci\ 



!£fc-ig-K-f^jfc_i)* 
\4\ 
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S2 











W\ 


■ \ci\ 


n 








\4\ 





[(o-l -CTj J ) + (o-i -crj2 ) + ■ ■ ■ + ((Tl -o-j^_ t 



(37) 



with 1 < ji < j2 < • • • < jk-i < 2A^. Thus, the time evolution of Di,ji,j2,- ;jk-i 
is mainly determined by the exponential law 



-^»Jl5j2viJfc — 1 



(38) 



Similarly, we deduce that the determinants that do not contain the first column 
of matrix C (!3T|) have the form 



D 



Jl 31 

\ci\ 141 



Rl Rl Ri 

^51 g-Cci-o-jj^)* ^j2.g-(cri-o-j2)< . . . ^fc g-(o'l-'^jfc)* 

Rl Rl Rl 



21^g-(0'l-0'jl)t f:^g-(0'l-0'j2)* ■ ■ ■ ^fc g-fgl-CTl,.)* 

I I J, I , 

"^1 



■^1 







cl 












f 




c2 














R" 



7^e 
1 



-[(cri-CTjJ + (cri-(Tj2)H h(o-i-crj^_ J + (cri-(Tjj.)]t 



(39) 



with 1 < J 1 < j2 < • • • < jfc-i < jfc < 2A^. Thus, the values of these determi- 
nants also tend to zero following an exponential law 



-[(cri-(Tjj) + (o-i-(Tj2)H l-(o-l-o-jj._j) + (o-l-crjj.)]t 



(40) 



Clearly, from all determinants appearing in the definition of Sk, flMl) . the one 
that decreases the slowest is the one containing the first k columns of matrix 
C in (EH): 

Di,2,3,...,k oc g-[K-<^2)+K-'^3)+-+K-'^fc)]* . (41) 



All other determinants appearing in equations (!38|) and fHOj) tend to zero faster 
than -Di,2,3,...,A; since the quantities in their exponentials are smaller or equal 
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to the exponent in (HTj) . We, therefore, conclude that the rate of decrease of 
Sk is dominated by (14T]) . yielding the approximation 

Sk{t) OC e-[('^l-'^2) + {ai-a3) + - + (ai-a,)]t _ (^42) 

Furthermore, since the norm f|T8l) of the /c-vector WiAw2A---AWkis expected 
to evolve in a similar way as 5*^, we conclude that GALI^ tends to zero in the 
same manner as above, i.e. 

GALIfe(t) OC e-[('^i-'^2)+('^i-'^^^)+-+('^i-'^*)l* . (43) 

We note here that in [17], where it was shown theoretically that SALI tends 
exponentially to zero for chaotic orbits as SALI(t) oc exp{ — (cti — o"2)t} (which 
is equivalent to equation fH3l) for k = 2), equation fl29|) was also retrieved, 
although it was wrongly assumed that the LCEs are related to the eigenvalues 
of matrix M of the variational equations f[T^ . 

In the previous analysis we assumed that o"i > cr2 so that the norm of each 
deviation vector can be well approximated by equation (130|) . If the first m 
Lyapunov exponents, with 1 < m < k, are equal, or very close to each other, 
i.e. 0"! ~ (72 — ■ • • — o"m equation (H3l) becomes 

GALIfc(t) OC e-[('^i-'^-+i)+('^i-'^™+2)+-+('^i-'^'=)l* , (44) 

which still describes an exponential decay. However, for k < m < N the GALIfe 
does not tend to zero as there exists at least one determinant of the matrix C 
that does not vanish. In this case, of course, one should increase the number 
of deviation vectors until an exponential decrease of GALIfc is achieved. The 
extreme situation that all cxj = corresponds to motion on quasiperiodic tori, 
where all orbits are regular and is described below. 



4-2 The evaluation of GALI for regular orbits 

As is well-known, regular orbits of an N degree of freedom Hamiltonian sys- 
tem flTUl) typically lie on A^-dimensional tori. If such tori are found around a 
stable periodic orbit, they can be accurately described by N formal integrals of 
motion in involution, so that the system would appear locally integrable. This 
means that we could perform a local transformation to action-angle variables, 
considering as actions Ji, J2, ■ ■ ■ , Jn the values of the formal integrals, so 
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that Hamilton's equations of motion, locally attain the form 



J^ = 

i = l,2,...,N. (45) 

9i = u!i{Ji, J2, . . . , Jn) 



These can be easily integrated to give 

J At) = Jin 

z = l,2,...,iV, (46) 

— ^io + ^i{Jio, J20, ■ ■ ■ , Jm) t 
where J^, 9io, i = 1,2, ... ,N are the initial conditions. 

By denoting as ^i, rji, i = 1,2, N small deviations of Jj and 9i respectively, 
the variational equations of system (jl5|) . describing the evolution of a deviation 
vector are 

6 = 

t = l,2,...,N, (47) 



where 

_ duji 



Jo 



z,j = l,2,...,N, (4J 



and Jo = (Jio, J20, ■ ■ ■ , Jm) = constant, represents the A^-dimensional vector 
of the initial actions. The solution of these equations is: 

r , t = l,2,...,N. (49) 



From equations (H9|) we see that an initial deviation vector w{0) with coordi- 
nates 6(0), 2 = 1, 2, . . . , in the action variables and ?7j(0), i = 1, 2, . . . , in 
the angles, i. e. w{0) = (6(0), 6(0), • • ., ^Ar(O), ?7i(0), 772(0), . . .,r]NiO)), evolves 
in time in such a way that its action coordinates remain constant, while its 
angle coordinates increase linearly in time. This behavior implies an almost 
linear increase of the norm of the deviation vector. To see this, let us assume 
that initially this vector td{0) has unit magnitude, i. e. 

N N 

E6(o)^ + E^^(o)' = i (50) 
1=1 1=1 
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whence the time evolution of its norm is given by 



Kt)ll= 1 + 



N / N 

EE ^^.0(0) 

z=l \j=l 



N I N 

2EU^(0)E^^.O(0) 

i=l \ j=l 




while the normalized deviation vector becomes 
1 



w{t) 



\w 



6(0),...,eiv(0),r7i(0) + 



E^i.O(O) 

3=1 



N 



E^^iO(O) 



t .(52) 



Since the norm (15T!) of a deviation vector, for t large enough, increases practi- 
cally linearly with t, the normalized deviation vector (l52l) tends to fall on the 
tangent space of the torus, since its coordinates perpendicular to the torus 
(i. e. the coordinates along the action directions) vanish following a rate. 
This behavior has already been shown numerically in the case of an integrable 
Hamiltonian of 2 degrees of freedom in [TCj. 



Using as a basis of the 2A^-dimensional tangent space of the Hamiltonian 
flow the 2N unit vectors {vi,V2, ■ ■ ■ ,V2n}, such that the first of them, 
Vi,V2, . . . ,vn correspond to the action variables and the remaining ones, 
■Oat+i, ■OAr+2, . . . , t'2Ar to the coujugatc angle variables, any unit deviation 



vector Wi, i = 1,2, . . . can be written as 



Wi{t) 



\w{t)\\ 



N 



N 



N 



E €(0) % + E v]io) + E ^kjQm 



k=l 



(53) 



We point out that the quantities 



i,j = 1,2. . .,N, in (jH]), depend only 
on the particular reference orbit and not on the choice of the deviation vector. 
We also note that the basis iti, i = 1,2, ... , 2N depends on the specific torus 
on which the motion occurs and is related to the usual vector basis Cj, z = 
1,2, ... , 2N of equation fll3p . through a non-singular transformation, similar 
to the one of equation fl35p . having the form: 



V1V2... V2N 



ei 62 



62N 



(54) 



with To denoting the transformation matrix. The basis {ei, 62, ... , e2iv} is used 
to describe the evolution of a deviation vector with respect to the original g^. 
Pi i = 1,2, ... ,N coordinates of the Hamiltonian system (fTOl) . while the basis 
{vi,V2, . . . ,V2n} is used to describe the same evolution, if we consider the 
original system in action-angle variables, so that the equations of motion are 
the ones given by (l45l) . 
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At this point we make the following remark: If the initial deviation vector 
already lies in the tangent space of the torus it will remain constant for all 
time! Indeed, taking for the initial conditions of this vector 

e.(0) = 0, 2 = 1,2,. ..,iV , (55) 

with 

N 

T.V^iOr = l, (56) 

i=l 

we conclude from equation fH9|) that 

, V^it)=V^iO)■ (57) 



i.e. the deviation vector remains unchanged having its norm always equal to 
1. In particular, such a vector has the form 

w{t) = (0, 0, . . . , 0, r/i(0), r/2(0), . . . , ^0)) . (58) 



Let us now study the case of k, general, linearly independent unit deviation 
vectors {wi,W2, ■ ■ ■ ,Wk} with 2 < k < 2N. Using as vector basis the set 
{vi,V2, . . .,V2n} we get: 



D 



ViV2 ... V2N 



(59) 



If no deviation vector is initially located in the tangent space of the torus, 
matrix D has the form 



D = [d. 



n 



k 

m=l 



em 

ef(o) 



a(0) r7K0) + E^=i^ime(0)t 
e^(0) r7?(0) + E^=i^i,ne(0)t 



r?^(0)+E™=1^7V,nd(0)t 

r?^(o)+E^):=i^7v™e(o)t 



er(o) ■ ■ ■ im vm + e:;=i ^i^e(o)t ■ ■ ■ r^m + ^N^imt 



where i = 1,2, . . . , k and j = 1,2, . . . , 2N. Recalling our earlier discussion (see 
fl50l) - fl53|) ). we note that the norm of vector Wi{t) for long times, grows linearly 
with t as 



Mi{t) = \\iUi{t)\\ (xt. 



(61) 
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Defining then by ^^''^ and rj'i the k x 1 column matrices 



0,k 



r/KO) ^.'(0) . . . ryf (0) 



(62) 



the matrix D of ( !60l) assumes the much simpler form 



D(t) 



^1 



k I x^N iO,k, 



1% + E»=l ^NiiTt 



0,k, 



nf=iM,(t) 



D°''-(^)- (63) 



Suppose now that we have m linearly independent deviation vectors, with 
m < k and m < N, initially located in the tangent space of the torus and let 
them be the first m deviation vectors in equation (159|) . This implies, in the 
above notation, that the vectors in (|63l) now have the form 



'.'m,k 



oo...oc'"+^(o)er^(o)...ef(o) 



m+2/ 



(64) 



where the first superscript, m, refers to the number of first components being 
equal to zero. Thus, the matrix D of (l63l) in this case reads 



D(t) 



m,k. 



i.m,k. 



D™''^(t), (65) 



nfrrM^+,(t) 



where the first superscript of D'"'^(t) in equations (l63l) and (165!) has an analo- 
gous meaning as in the We note that for = m we define HLi Mm+i{t) = 
1. 



Using again equation ( lA.lSp . we write the wedge product of the k normalized 
deviation vectors as 



wi{t) Aw2{t) A- ■ ■ AWk{t) = 

l<n<i2<--<«fe<2Af 



d2ii d2i2 ■ ■ ■ d2i^ 



dkil dki2 ' ' ' dj^i^ 



Ui^AUi^A--- A tiife.(66) 
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and introduce the analogous quantity 



Sh. 



E 

l<n<i2<-<ife<2Af 



diii dii^ ■ ■ ■ dii^, 
d2ii d2i2 ■ ■ ■ d2if. 



dkil d}^i2 ■ ■ ■ C^fcij. 



1/2 



(67) 



as in the case of chaotic orbits, see fl33l) and flM|) respectively. 

As we have already explained, the k deviation vectors will eventually fall on 
the iV-dimensional tangent space of the torus on which the motion occurs. Of 
course, if some of them are already located in the tangent space, at t = 0, they 
will remain there forever. In their final state, the deviation vectors will have 
coordinates only in the A^-dimensional space spanned by vn+i, vn+2, ■ ■ ■ , V2n- 
Now, if we start with 2 < k < N general deviation vectors there is no particu- 
lar reason for them to become linearly dependent and their wedge product will 
be different from zero, yielding 5*^ and GALI^ which are not zero. However, 
if we start with N < k < 2N deviation vectors, some of them will necessarily 
become linearly dependent. Thus, in this case, their wedge product (as well 
as S'l^ and GALIfc) will be zero. 



We, therefore, need to examine in more detail the behavior of these 5*^. Since, 
in general, we choose the initial deviation vectors randomly (insisting only that 
they be linearly independent), the most common situation is that none of the 
initial deviation vectors is tangent to the torus. However, as we are not certain 
that this will always hold, let us suppose that < m < of our deviation 
vectors are initially in the tangent space of the torus. For 2 < k < N, this 
will make no difference, as the GALI^ tends to a non-zero constant. However, 
for N < k < 2N, GALI^ goes to zero by a power law and the fact that m 
vectors are already in the tangent space, at t = 0, may significantly affect the 
decay rate of the index. Thus, in such cases, the behavior of GALI needs to 
be treated separately. 



4-2.1 The case of m = tangent initial deviation vectors 

Let us consider first the most general case that no deviation vector is ini- 
tially tangent to the torus. In this case, the matrix D, whose elements appear 
in the definition of S*^, has the form given in equation (l63l) . So all determi- 
nants appearing in the definition of S"^ have as a common factor the quantity 
V riiLi ^j('^)) which, due to (ETj) . decreases to zero according to the power 
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law 



oc 



niiM,(t) t 



(6J 



In order to determine the precise time evolution of S'j., we search for the fastest 
increasing determinants of all the possible k x k minors of the matrix D*^'*^, in 
f lB51) . as time t grows. 

Let us start with k being less than or equal to the dimension of the tangent 
space of the torus, i. e. 2 < A; < A^. The fastest increasing determinants in 
this case are the N\/{k\{N — k)\) determinants, whose k columns are chosen 
among the last columns of matrix D*^''^: 



jl,j2,--;jk 



with 1 < ji < j2 < • • • < Jfc ^ N. Using standard properties of determinants, 
we easily see that the time evolution of is mainly determined by the 

behavior of determinants of the form 



^jimi^m^t UJ. 



Ijj 



tO,k tO,k . . . t 



0,fc 
rrife 



oct^(7 



where rrii G {1,2,..., A^}, i = 1,2, ... ,k, with ^ rrij, for all i ^ j. Thus, 
from (168|) and (ITOll we conclude that the contribution to the behavior of S"^ of 
the determinants related to ^^f,j2,-,jk ^'^ provide constant terms in (1671) . All 
other determinants appearing in the definition of 5^, not being of the form 
of A'?''^,- , contain at least one column from the first A^ columns of matrix 
D^'^^' and introduce in fl67|) terms that grow at a rate slower than t^, which 
will ultimately have no bearing on the behavior of GALIfe(t). To see this, let 
us consider a particular determinant of this kind 



0,A: 



SI ■ ■ ■ Vl+ l^i=l ^li^i ^ ■ ■ ■ Vk-m + Z-i=l ^^k-mi^i t 



(71) 



containing the first m columns of matrix D^''^, which are related to the action 
coordinates of the system, and the first k — m columns of the angle related 
columns of D^''^, with 1 < m < k. The first m columns of A^'^ are time 
independent. Using repeatedly a standard property of determinants, we easily 
see that the time evolution of is mainly determined by the time evolution 
of determinants of the form: 



?1 ?2 ^lil^h ^ ^2i2^i2 I ■ ■ ■ ^k-mik-rr^^ik-J- 



OC r-™, (72) 
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with ij G {m + 1, m + 2, . . . , A^}, j = 1,2, . . . , k — m and ij ^ ii, for all j ^ I. 
Thus, the contribution to the behavior of S"^ of determinants similar to A^*^ 
are terms proportional to t^'"^ /t^ = l/t"^ {1 < m < k), tending to zero as t 
grows. Since the k x k determinants appearing in the definition of S*^ involve 
both terms of the form fl69l) . growing as t'' and of the form fl7T|) . growing as 
fk-m^ the overall behavior of S*^ will be defined by determinants growing as 
t^, which when combined with f l68|) yields the important result 

GALIfc(t) ^ constant for 2 < k < N. (73) 



Next, let us now turn to the case of k deviation vectors with N < k < 2N. The 
fastest growing determinants are again those containing the last columns 
of the matrix D°'^: 



. 0,A; 

il J2,--- j'fe-JV 



1,2,. ..,Ar 



fO,fe 



N 



(74) 



with 1 < ji < j2 < • • • < jk~N < N. The first k-N columns of ^°[%^,„j^_^^i^2,...,n 
are chosen among the first columns of D'^''^ which are time independent. 
So there exist N\/{{k — N)\{2N — k)\) determinants of the form (1741) . which 
can be written as a sum of simpler k x k determinants, each containing in the 
position of its last columns r]^, i = 1,2, . . . ,N and/or columns of the form 
Uji^i'^t with i,j = 1,2, . . . , N. We exclude the ones where z = 1, 2, . . . , A^ 
appear more than once, since in that case the corresponding determinant is 
zero. Among the remaining determinants, the fastest increasing ones are those 
containing as many columns proportional to t as possible. 

^0,k 



Since t is always multiplied by the ^j' , and such columns occupy the first 

0,fc 

ii j2,---,ifc-iv,i 



k — N columns of A^''^- • i 2 ... at, ^ appears at most N — {k — N) = 2N — k 



times. Otherwise the determinant would contain the same ^^''^ column at least 
twice and would be equal to zero. The remaining k — {2N — k) — {k—N) = k — N 
columns are filled by the rjf each of which appears at most once. Thus, the 
time evolution of ^'jfj2,-,jk-NA,'2.-,N is mainly determined by determinants of 
the form: 

^0,fc ^O.fc 0,fc 0,fc tO,k. t^,k , ^ ^2N-k (yr\ 



with ii G {1, 2, . . . , A^}, / = 1,2,...,A^, ii 7^ ij, for all / 7^ j and m; G 
{1,2, ... ,A^}, / = 1,2, . . . ,2A^ - fc, mi ^ {ji,j2, • • • , jfc-7v}, mi 7^ mj, for 
all I 7^ j. So determinants of the form (17^ contribute to the time evolu- 
tion of by introducing terms proportional to = l/t^*^'^~^^. All 
other determinants appearing in the definition of S*^, not having the form of 
A^''^. 1 o AT, introduce terms that tend to zero faster than l/t^(^~^) 
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since they contain more than k — N time independent columns of the form 
^^'^^j i = 1,2, . . . , N. Thus S'j. and consequently GALI^ tend to zero following 
a power law of the form: 

GALIfc(t) cx --^ for iV < A; < 2N. (76) 



4-2.2 The case of m > tangent initial deviation vectors 

Finally, let us consider the behavior of GALI^ for the special case where m 
initial deviation vectors, with m < k and m < N, are located in the tangent 
space of the torus. In this case, matrix D, whose elements appear in the 
definition of 5*^, has the form given by (1651) . Thus, all determinants appearing 
in the definition of S'f. have as a common factor the quantity 1/ nfj™ Mm+i{t), 
which decreases to zero following a power law 



Proceeding in exactly the same manner as in the m = case above, we deduce 
that, in the case oi2 <k < N the fastest growing k^k determinants resulting 
from the matrix D"^''^ are of the form: 



UJ. 



'^k'^k — m ' 



t 



with ii G {1, 2, . . . , A^}, I = l,2,...,k with ii ^ ij for / ^ j, and ni G 
{1, 2, . . . , A^}, I = 1,2, . . . , k — m with ni ^ nj, for / ^ j. Hence, we conclude 
that the behavior of S'^., and consequently of GALI^ is defined by the behavior 
of determinants having the form of fl78|) which, when combined with flTTl) 
implies that 

GALIfc(t) ^ constant for 2 < k < N. (79) 



The case of A^ < A; < 2A^ deviation vectors, however, with m > initially 
tangent vectors, yields a considerably different result. Following entirely anal- 
ogous arguments as in the m = case, we find that, if m < A; — A^, S*^ 
and GALIfc evolve proportionally to f^-k^^k-m ^ ^^^2(k-N)-m_ ^^^lei 
hand, if m > k — N, one can show that the fastest growing determinant is 
proportional to t^"*". In this case, 5*^ and GALI^ evolve in time following a 
quite different power law: t'^-''^ /tk-m _ i^fk~N_ 

Summarizing the results of this section, we see that GALI^ for regular motion 
remains essentially constant when k < N, while it tends to zero for k > N 
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following a power law which depends on the number m {m < N and m < k) of 
deviation vectors initially tangent to the torus. In conclusion, we have shown 
that: 



constant if 2 < A; < 



GALIfc(t) oc < 



1 



if < A; < 2A^ and < m < A; - iV 



f.2(k-N)-m 

^ if < A; < 2A^ and m > A; - 



(80) 



k-N 



5 Numerical verification and applications 



In order to apply the GALI method to Hamiltonian systems and verify the 
theoretically predicted behavior of the previous sections, we shall use two 
simple examples with 2 (2D) and 3 (3D) degrees of freedom: the well-known 
2D Henon-Heiles system [18], described by the Hamiltonian 

H2 = lipl + pD + \{^^ + y^) + ^^y - \y^ (si) 



and the 3D Hamiltonian system: 



3 ^. 

^3 = 51 -^(^^i + Pi) + ^1^2 + g?g3, 
i=i ^ 



studied in 
energies H2 



9l[5] . We keep the parameters of the two systems fixed at the 



0.125 and = 0.09, with lui = 1, uj2 = and 



UJ3 



Vs. In 



order to illustrate the behavior of GALI^, for different values of A;, we shall 
consider some representative cases of chaotic and regular orbits of the two 
systems. 

Additionally, we shall study the higher-dimensional example of a 15D Hamil- 
tonian, describing a chain of 15 particles with quadratic and quartic nearest 
neighbor interaction, known as the famous Fermi-Pasta-Ulam (FPU) model 
150] 



1 



15 



15 



^15 = ^Ep.' + E 



=1 



where is the displacement of the ith particle from its equilibrium point and 
Pi is the conjugate momentum. This is a model we have recently analyzed in 
[39] and we shall use here the same values of the energy i^is = 26.68777 and 
/5 = 1.04 as in that study. 
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Fig. 1. (a) The evolution of Li(t) (solid curve), L2{t) (dashed curve) and Li{t)—L2{t) 
(dotted curve) for a chaotic orbit with initial conditions x = 0, y = —0.25, px = 0.42, 
Py = Qoi the 2D system (b) The evolution of GALI2, GALI3 and GALI4 of 

the same orbit. The plotted lines correspond to functions proportional to e"'^^* 
(solid hue), e^^'^i* (dashed line) and e"^'^!* (dotted hue) for cji = 0.047. Note that 
the t-axis is linear. The evolution of the norm of the deviation vector 'w{t) (with 
||tt;(0)|| = 1) used for the computation of Li{t), is also plotted in (b) (gray curve). 

5.1 A 2D Hamiltonian system 



Let us consider first a chaotic orbit of the 2D Hamiltonian (IHTI) . with initial 
conditions x = 0, ?/ = —0.25, = 0.42, py = 0. In figure [T](a) we see the time 
evolution of Li{t) of this orbit. The computation is carried out until Li{t) 
stops having large fluctuations and approaches a positive value (indicating 
the chaotic nature of the orbit), which could be considered as a good approx- 
imation of the maximal LCE, ai. Actually, for t ^ 10^, we find cti ^ 0.047. 

We recall that 2D Hamiltonian systems have only one positive LCE cxi, since 
the second largest is (T2 = 0. It also holds that cxs = —(T2 and = — cti and 
thus formula (H3|) . which describes the time evolution of GALI^ for chaotic 
orbits, gives 

GALl2(t) oc e""^*, GALl3(t) oc e^^"^*, GALl4(t) oc e"^"^*. (84) 



In figure [T](b) we plot GALI^, A; = 2, 3, 4 for the same chaotic orbit as a func- 
tion of time t. We plot t in linear scale so that, if (18^ is valid, the slope of 
GALI2, GALI3 and GALI4 should approximately be — ai/lnlO, — 2(Ti/lnlO 
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and —4(71 /In 10 respectively. From figure [T](b) we see that lines having pre- 
cisely these slopes, for ai = 0.047, approximate quite accurately the computed 
values of the GALIs. The biggest deviation between the theoretical curve and 
numerical data appears in the case of GALI4 where the theoretical prediction 
underestimates the decaying rate of the index, but even in this case the differ- 
ence does not appear too significant. Note, however, the important difference 
in the times it takes to decide about the chaotic nature of the orbit: Waiting 
for the maximal LCE to converge in figure [U^a), one needs more than 10^ t ime 
units, while, as we see in figure [T](b), the GALI^'s provide this information in 
less than 400 time units! 

We also note that, plotting in this example the evolution of the quantity 
(with ||w(0)|| = 1), which is used to determine Li(t) in ([1]) and 
is practically identified with the Fast Lyapunov Indicator (FLI), we obtain 
in figure [U^b) a graph similar to that of GALl2(t). This is not surprising, 
as both and GALl2(t) tend exponentially to zero following a decay 

proportional to e"'^^* (see equations fl30l) and fl8^ ). From the results of figure 
Wip) we see that the different plotted quantities reach the limit of computer's 
accuracy (10~^^) at different times and in particular GALI2 a.tt^ 800, GALI3 
at t ^ 400, GALI4 at t ^ 150 and ||w(t)||-i at t ^ 720. The CPU time needed 
for computing the evolution of the indices up to these times were: 0.220 sec 
for \\w{t)\\-\ 0.295 sec for GALI2, 0.165 sec for GALI3 and 0.070 sec for 
GALI4 respectively. Thus, in this case also, it is clear that the higher order 
GALIfc (with k > 2) can identify the chaotic nature of an orbit faster than the 
methods of the maximal LCE, the FLI or the SALI (equivalent to GALI2, see 
below) . 

It is interesting to remark at this point (as mentioned in section 14. ip . that 
the accuracy of the exponential laws ([HID is due to the fact that the local 
Lyapunov exponents cease to fluctuate significantly about their limit values, 
after a relatively short time interval. To see this, we have plotted in figure 
[U^a), the two nonnegative local Lyapunov exponents Li(t), L2(t), as well as 
their difference. Note that Li{t) — L2{t) begins to be well approximated by 
cTi — (T2 = 0"! already for times t of order 10^ units. A similar behavior of 
such Li(t) — Li(t), i = 2,3, . . . ,2N differences are observed for the other 
Hamiltonians we studied in this paper having 3 or more degrees of freedom. 

As explained in detail in Appendix [Bl GALI2 practically coincides with SALI 
in the case of chaotic orbits. This becomes evident from figure |2] where we plot 
the absolute difference between GALI2 and SALI of the chaotic orbit of figure 
[T] as a function of time t. The two indices practically coincide after about 
t ~ 300 units, since their difference is at the limit of computer's accuracy 
(10~^^), although their actual values are of order 10~^ (see figure [Tt^b)). 

Let us now study the behavior of GALI^ for a regular orbit of the 2D Hamil- 
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Fig. 2. The absolute difference between GALI2 and SALI of the chaotic orbit of 
figure [Das a function of time t. 

tonian (IHTl) . From (IHOl) it follows that in the case of a Hamiltonian system 
with N = 2 degrees of freedom GALI2 will always remain different from zero, 
while GALI3 and GALI4 should decay to zero following a power law, whose 
exponent depends on the number m of deviation vectors that are initially tan- 
gent to the torus on which the orbit lies. Now, for a regular orbit of the 2D 
Hamiltonian ( IHTj) and a random choice of initial deviation vectors, we expect 
the GALI indices to behave as 

GALl2(t) oc constant, GALl3(t) oc \, GALl4(t) oc ^. (85) 



A simple qualitative way of studying the dynamics of a Hamiltonian system is 
by plotting the successive intersections of the orbits with a Poincare Surface 
of Section (PSS) [45j. In 2D Hamiltonians, the PSS is a two dimensional plane 
and the points of a regular orbit (which lie on a torus) fall on a smooth closed 
curve. This property allows us to choose initial deviation vectors tangent to a 
torus in the case of system (IHTI) . In particular, we consider the regular orbit 
with initial conditions x = 0, y = 0, = 0.5, Py = 0. In figure [31 we plot the 
intersection points of this orbit with the PSS defined by x = (panel (a)) and 
y = (panel (b)). From the morphology of the two closed curves of figure [3l 
it is easily seen that deviation vectors ei = (1, 0, 0, 0) and 64 = (0, 0, 0, 1) are 
tangent to the torus. 

In Figure HI we plot the time evolution of SALI, GALI2, GALI3 and GALI4 
for the regular orbit of figure [31 for various choices of initial deviation vectors. 
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Fig. 3. The Poincare Surface of Section (PSS) defined by (a) x = and (b) y = 
of the regular orbit with initial conditions x = 0, y = 0, Px = 0.5, py = for the 
Henon-Heiles system (jSip . 

In figure Hl^a) the initial deviation vectors are randomly chosen so that none of 
them is tangent to the torus. In this case SALI and GALI2 fluctuate around 
non-zero values, while GALI3 and GALI4 tend to zero following the theo- 
retically predicted power laws, see (!85l) . In figure |l](b) we present results for 
the indices when we have m = 1 initial deviation vector tangent to the torus 
(in particular vector ei). In this case the indices evolve as predicted by (IHOl) . 
i. e. SALI and GALI2 remain practically constant, while GALI3 oc 1/t and 
GALI4 oc Finally, in figure |l](c) we have plotted our results using m = 2 
initial deviation vectors tangent to the torus (vectors ei and 64). Again the 
predictions of ( !80l) are seen to be valid since GALI3 oc l/t and GALI4 oc 

The different behavior of SALI (or GALI2) for regular and chaotic orbits has 
already been successfully used for discriminating between regions of order and 
chaos in various dynamical systems [17 36 4 d|41|l42|l43l|44] . For example, by 
integrating orbits whose initial conditions lie on a grid, and by attributing 
to each grid point a color according to the value of SALI at the end of a 
given integration time, one can obtain clear and informative pictures of the 
dynamics in the full phase space of several Hamiltonian systems of physical 
significance [T7l36|43j . 

Figures [U^b) and H] clearly illustrate that GALI3 and GALI4 tend to zero both 
for regular and chaotic orbits, but with very different time rates. We may use 
this difference to distinguish between chaotic and regular motion following a 
different approach than SALI or GALI2. Let us illustrate this by considering 
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Fig. 4. Time evolution of SALI (gray curves), GALI2, GALI3 and GALI4 for the 
regular orbit of figure [3] in log-log scale for different values of the number m of 
deviation vectors initially tangent to the torus: (a) m = 0, (b) m = 1 and (c) m = 2. 
We note that in panel (a) the curves of SALI and GALI2 are very close to each 
other and thus cannot be distinguished. In every panel, dashed lines corresponding 
to particular power laws are also plotted. 



the computation of GALI4: From fl84l) and fl85l) . we expect GALI4 oc e~''°"i* 
for chaotic orbits and GALI4 oc for regular ones. These time rates imply 
that, in general, the time needed for the index to become zero is much larger 
for regular orbits. Thus, instead of simply registering the value of the index at 
the end of a given time interval (as we do with SALI or GALI2), let us record 
the time, tj^, needed for GALI4 to reach a very small threshold, e. g. 10~^^, 
and color each grid point according to the value of tth- 
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Fig. 5. Regions of different values of the time tth needed for GALI4 to become less 
than 10^^^ on the PSS defined by x = of the 2D Henon-Heiles Hamiltonian (I8ip . 

The outcome of this procedure for the 2D Henon-Heiles system ( IHTi) is pre- 
sented in figure [5l Each orbit is integrated up to t = 500 units and if the value 
of GALI4, at the end of the integration is larger than 10~^^ the corresponding 
grid point is colored by the light gray color used for tth > 400. Thus we can 
clearly distinguish in this figure among various 'degrees' of chaotic behavior in 
regions colored black or dark gray - corresponding to small values of tth ^ and 
regions of regular motion colored light gray, corresponding to large values of 
tth- At the border between them we find points having intermediate values of 
tth which belong to the so-called 'sticky' chaotic regions. Thus, this approach 
yields a very detailed chart of the dynamics, where even tiny islands of sta- 
bility can be identified inside the large chaotic sea. We note that for every 
initial condition the same set of initial deviation vectors was used, ensuring 
the same initial value of GALI4 for all orbits and justifying the dynamical 
interpretation of the color scale of figure [51 



5.2 A 3D Hamiltonian system 

Let us now study the behavior of the GALIs in the case of the 3D Hamiltonian 
(1521) . Following [ISfS] the initial conditions of the orbits of this system are 
defined by assigning arbitrary values to the positions gi, q2, qs, as well as the 
so-called 'harmonic energies' Ei, E2, related to the momenta through 

, ^ 

P^ = \ — , ^ = 1,2,3. 86 

V 
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Fig. 6. (a) The evolution of Li{t), L2{t) for the chaotic orbit with initial condition 
= 52 = 93 = 0, Ei = E2 = E3 = 0.03 of the 3D system ([82]). (b) The evolution of 
GALIfe with A; = 2, . . . , 6 of the same orbit. The plotted lines correspond to functions 
proportional to e'^'^i-'^^)*^ g-(2ai-<72)t^ g-(3ai-<72)i^ g-4<7it ^-Qa,t £qj. ^ Q 

(72 = 0.008. Note that the t-axis is linear. 

Chaotic orbits of 3D Hamiltonian systems generally have two positive Lya- 
punov exponents, ai and (72, while a-^ = 0. So, for approximating the behavior 
of GALIs according to (l43l) . both ai and a2 are needed. In particular, (H3l) 
gives 

GALl2(t) oc e-("i-"^)*, GALl3(t) oc e-^^-i— 2)*, GALl4(t) oc g-^^-i— 2)*^ 
GALl5(t) oc e-^"i*, GALl6(t) oc e"^"!*. 



Let us consider the chaotic orbit with initial conditions qi = q2 = = 
0, El = E2 = E3 = 0.03 of the 3D system (JHSl). We compute ai, ^2 for 
this orbit as the long time limits of the Lyapunov exponent quantities Li{t), 
L2(t), applying the technique proposed by Benettin et al. [5]. The results are 
presented in figure [6](a). The computation is carried out until Li{t) and L2{t) 
stop having large fluctuations and approach some positive values (since the 
orbit is chaotic), which could be considered as good approximations of their 
limits (Ji, (72. Actually for t ~ 10^ we have (Xi ^ 0.03 and (72 ~ 0.008. Using 
these values as good approximations of cti, (72 we see in figure [6]^b) that the 
slopes of all GALIs are well reproduced by fl87|) . 

Next, we consider the case of regular orbits in our 3D Hamiltonian system. 
In the general case, where no initial deviation vector is tangent to the torus 
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Fig. 7. (a) The evolution of Li(t) for the regular orbit with initial condition 
= 52 = 93 = 0, El = 0.005, E2 = 0.085, = of the 3D system ([82]). (b) 
The evolution of GALI^ with k = 2, ... ,6 of the same orbit. The plotted lines 
correspond to functions proportional to p-, ^ and ^. 

where the regular orbit lies, the GALIs should behave as: 

GALl2(t) oc constant, GALl3(t) oc constant, GALl4(t) oc ^, 
GALl5(t) oc ^, GALl6(t) oc ^. 



according to (IHOl) . In order to verify expression (188|) we shall follow a specific 
regular orbit of the 3D system (!82|) with initial conditions qi = q2 = Qs = 0, 
El = 0.005, E2 = 0.085, E3 = 0. The regular nature of the orbit is revealed 
by the slow convergence of its Li(t) to zero, implying that ai = 0, see figure 
[Tl^a). In figure [7](b), we plot the values of all GALIs of this orbit with respect 
to time t. From these results we see that the different behaviors of GALIs are 
very well approximated by formula fl88l) . 

From the results of figures M and U\ therefore, we conclude that in the case of 
3D Hamiltonian systems not only GALI2, but also GALI3 has different behav- 
ior for regular and chaotic orbits. In particular GALI3 tends exponentially to 
zero for chaotic orbits (even faster than GALI2 or SALI), while it fluctuates 
around non-zero values for regular orbits. Hence, the natural question arises 
whether GALI3 can be used instead of SALI for the faster detection of chaotic 
and regular motion in 3D Hamiltonians and, by extension, whether GALI^, 
with k > 3, should be preferred for systems with > 3. The obvious compu- 
tational drawback, of course, is that the evaluation of GALI^ requires that we 
numerically follow the evolution of more than 2 deviation vectors. 
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First of all, let us point out that the computation of SALT, applying (Q, 
is slightly faster than GALI2, for which one needs to evaluate several 2x2 
determinants. For example, for orbits of the 3D Hamiltonian fl82l) the CPU 
time needed for the computation of SALI for a fixed time interval t, was about 
97% of the CPU time needed for the computation of GALI2 for the same time 
interval. Although this difference in not significant, we prefer to compute SALI 
instead of GALI2 and compare its efficiency with the computation of GALI3. 

It is obvious that the computation of GALI3 for a given time interval t needs 
more CPU time than SALI, since we follow the evolution of three deviation 
vectors instead of two. This is particularly true for regular orbits as the index 
does not become zero and its evolution has to be followed for the whole pre- 
scribed time interval. In the case of chaotic orbits, however, the situation is 
different. Let us consider, for example, the chaotic orbit of figured The usual 
technique to characterize an orbit as chaotic is to check, after some time inter- 
val, if its SALI has become less than a very small threshold value, e. g. 10~^. 
For this particular orbit, this threshold value was reached for t ~ 760. Adopt- 
ing the same threshold to characterize an orbit as chaotic, we find that GALI3 
becomes less than 10~^ after t ~ 335, requiring only as much as 65% of the 
CPU time needed for SALI to reach the same threshold! 

So, using GALI3 instead of SALI, we gain considerably in CPU time for chaotic 
orbits, while we lose for regular orbits. Thus, the efficiency of using GALI3 
for discriminating between chaos and order in a 3D system depends on the 
percentage of phase space occupied by chaotic orbits (if all orbits are regu- 
lar GALI3 requires more CPU time than SALI). More crucially, however, it 
depends on the choice of the final time, up to which each orbit is integrated. 
As an example, let us integrate, up to t = 1000 time units, all orbits whose 
initial conditions lie on a dense grid in the subspace ^3 = P3 = 0, p2 > of a 
4-dimensional PSS, with gi = of the 3D system flS^ . attributing to each grid 
point a color according to the value of GALI3 at the end of the integration. 
If GALI3 of an orbit becomes less than 10~® for t < 1000 the evolution of the 
orbit is stopped, its GALI3 value is registered and the orbit is characterized 
as chaotic. The outcome of this experiment is presented in figure [HI 

We find that 77% of the orbits of figure [Dare characterized as chaotic, having 
GALI3 < 10~^. In order to have the same percentage of orbits identified as 
chaotic using SALI (i. e. having SALI < 10~^) the same experiment has to 
be carried out for t = 2000 units, requiring 53% more CPU time. Due to 
the high percentage of chaotic orbits, in this case, even when the SALI is 
computed for t = 1000 the corresponding CPU time is 12% higher than the 
one needed for the computation of figure [H while only 55% of the orbits 
are identified as chaotic. Thus it becomes evident that a carefully designed 
application of GALI3 - or GALI^ for that matter - can significantly diminish 
the computational time needed for a reliable discrimination between regions 
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Fig. 8. Regions of different values of the GALI3 on the subspace qs = Ps = 0, P2 ^ 
of the 4-dimensional PSS gi = of the 3D system ([82]) at t = 1000. 

of order and chaos in Hamiltonian systems with N > 2 degrees of freedom. 



5.3 A multi-dimensional Hamiltonian system 



Let us finally turn to a much higher-dimensional Hamiltonian system hav- 
ing 15 degrees of freedom, i. e. the one shown in (183|) . With fixed boundary 
conditions 

goW = gieW = 0, Vt, (89) 



it is known that there exists, for all energies, iJis = E, a simple periodic orbit, 
satisfying [?T|5^ 

q2i{t) = 0, q2i-i{t) = -q2i+i{t) = q{t), 2 = 1,2,..., 7, (90) 



where q{t) = q{t + T) obeys a simple nonlinear equation admitting Jacobi 
elliptic function solutions. For the parameter values Hi^ = 26.68777 and /? = 
1.04 used in an earlier study [39], we know that this orbit is unstable and has 
a sizable chaotic region around it. As initial conditions for (!90l) we take 

q{Q) = 1.322 and pi{0) = 0, i = 1, 2, . . . , 15. (91) 



First, we consider a chaotic orbit which is located close to this periodic so- 
lution, by taking as initial conditions qi{0) = q{0), ^3(0) = 57(0) = gii(O) = 
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Fig. 9. (a) The evolution of Li{t), L2{t), L3{t) and L4(t) for a chaotic orbit of the 
15D system ([83]). (b) The evolution of GALI2, GALI3 and GALI4 for the same orbit. 
The plotted lines correspond to functions proportional to e~^°"^~°'^^^ , ^-C^^^i- (^2-0-^)1 
and e-(3'^i-'^2-<T3-^4)i^ fo^ ai = 0.132, era = 0.117, 0-3 = 0.104, 0-4 = 0.093. Note 
that the t-axis is linear. 

-q{0) + 10-\ g5(0) = g9(0) = ^15(0) = g(0) - 10"^ gai = for z = 1, 2, . . . , 7 
and Pi{0) = for z = 1,2, ...,14, pi5(0) = 0.00323. The chaotic nature of 
this orbit is revealed by the fact that its maximal LCE is positive (see figure 
[9]^a)). In fact, from the results of figure [9](a) we deduce reliable estimates of the 
system's four largest Lyapunov exponents: cti ~ 0.132, (T2 ~ 0.117, ~ 0.104 
and CT4 ~ 0.093. Thus, we have a case where several LCEs have positive values, 
the largest two of them being very close to each other. The behavior of the 
GALIs is again quite accurately approximated by the theoretically predicted 
exponential laws (H3l) . This becomes evident by the results presented in figure 
Mjo), where we plot the time evolution of GALI2, GALI3 and GALI4 as well as 
the exponential laws that theoretically describe the evolution of these indices. 
In this case, GALI2 does decay to zero relatively slowly since ui and (T2 have 
similar values and hence, using GALI3, GALI4 or a GALI of higher order, one 
can determine the chaotic nature of the orbit much faster. 

It is worth mentioning that P3|) describes much more accurately the evolution 
of GALIfe when the orbit we wish to study is very close to the unstable periodic 
solution fl90|) itself. This is due to the fact that in that case, the LCEs are 
directly related to the eigenvalues of the monodromy matrix associated with 
the variational equations of this unstable periodic orbit, see equation (^^. 
In fact, for our choice of parameters, this matrix has two equal pairs of real 
eigenvalues with magnitude greater than one, while all other eigenvalues lie 
on the unit circle in the complex plane. As a consequence, the orbit has two 
nearly identical positive Lyapunov exponents (as well as their two negative 
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Fig. 10. (a) The evolution of Li{t), L2{t), L^{t) and L4^{t) for an orbit which is very 
close to the unstable periodic orbit ([9T]) of the 15D system (f83]l . (b) The evolu- 
tion of GALI2, GALI3 and GALI4 of the same orbit. The plotted lines correspond 
to functions proportional to e'^''^-''^^^ ^-{2^7^-02-0^)1 ^-(301-02-03-0^)1^ 

cji = 0.3885, (T2 = 0.3883, (T3 = 0, CJ4 = 0. Note that the t-axis is hnear. 

counterparts), while all other exponents are zero. This is shown in figure [TOlfa). 
where we plot the evolution of the Li{t) for i = 1,2,3,4, whose limits for 
t — s> 00 are the 4 largest Lyapunov exponents. From these results we deduce 
(Ti ~ 0.3885, (J2 ~ 0.3883, while the decrease of L^it) and L^it) to zero 
indicate that o"3 = (74 = 0. In figure [TOlfb) we now observe that GALI2 remains 
practically constant for this particular time interval (actually it decreases to 
zero extremely slowly following the exponential law e"*^'^^"'^^^* = e~° '"^°^*). On 
the other hand, GALI3 and GALI4 decay exponentially to zero following the 
laws, GALI3 oc e-(2'^i-'^2-fT3)t^ GALI3 oc e~^^''^-''^~''^-''^^\ given by equation 



6 Discussion and conclusions 



In this paper we have introduced and applied the Generalized Alignment In- 
dices of order k (GALI^) as a tool for studying local and global dynamics in 
conservative dynamical systems, such as Hamiltonian systems of degrees 
of freedom, or 2A^-dimensional symplectic maps. We have shown that these 
indices can be successfully employed not only to distinguish individual orbits 
as chaotic or regular, but also to efficiently chart large domains of phase space, 
characterizing the dynamics in the various regions by different behaviors of the 
indices ranging from regular (GALI^s are constant or decay by well-defined 
power laws) to chaotic (GALI^s exponentially go to zero). 
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A different approach than simply calculating the maximal Lyapunov exponent 
is to compute the so-called Smaller Alignment Index (SALT), following the 
evolution of two initially different deviation vectors. This approach has been 
used by several authors and has proved quite successful, as it can determine 
the nature of the dynamics more rapidly, reliably and efficiently than the 
maximal LCE. In the present paper, motivated by the observation that the 
SALI is in fact proportional to the 'area' of a parallelogram, having as edges 
the two normalized deviation vectors, we have generalized SALI by defining a 
quantity called GALI^, representing the 'volume' of a parallelepiped having as 
edges k > 2 initially linearly independent unit deviation vectors. In practice, 
GALIfc is computed as the 'norm' of the 'exterior' or wedge product of the k 
normalized deviation vectors. 

For the numerical evaluation of GALI^, we need to compute the reference 
orbit we are interested in from the fully nonlinear equations of the system, as 
well as follow the time evolution of k deviation vectors, solving the (linear) 
variational equations about the orbit. How many such vectors should we take? 
Since the phase space of the dynamical system is 2A^-dimensional, k should 
be less than or equal to 2N, otherwise GALIfc will be equal to zero already 
from the start. However, even though we may choose our deviation vectors 
initially linearly independent, they may become dependent as time evolves, in 
which case the phase space 'volume' represented by GALI^ will vanish! This 
is precisely what happens for all > 2 if our reference orbit is chaotic, and 
also if it is regular and k > N, but at very different time rates. 

In particular, we showed analytically and verified numerically in a number 
of examples of Hamiltonian systems that for chaotic orbits GALI^ tends ex- 
ponentially to zero following a rate which depends on the values of several 
Lyapunov exponents (see equation fH5]) ). On the other hand, in the case of 
regular orbits, GALI^ with 2 < k < N fiuctuates around non-zero values, 
while, for N < k < 2N, it tends to zero following a power law (see equation 
(IHO!)). The exponent of the power law depends on the values of k and A^, as well 
as on the number m of deviation vectors that may have been chosen initially 
tangent to the torus on which the orbit lies. 

Clearly, these different behaviors of the GALI^, can be exploited for the rapid 
and accurate determination of the chaotic versus regular nature of a given 
orbit, or of an ensemble of orbits. Varying the number of deviation vectors 
(and bringing more LCEs into play), we can, in fact, achieve high rates of 
identification of chaotic regions, in a computationally advantageous way. Sec- 
ondly, regular motion can be identified by the index being nearly constant for 
small k, while, when k exceeds the dimension of the orbits' subspace, GALI^ 
decays by well-defined power laws. This may help us identify, for example, 
cases where the motion occurs on cantori of dimension d < N (see e.g. |45j ) 
and the orbits become 'sticky' on island chains, before turning truly chaotic 
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and exponential decay takes over. 

We have also studied on specific Hamiltonians with N > 2 the computational 
efficiency of the GALI^. One might suspect, of course, that the best choice 
would be GALIiv since this is the index that exhibits the most different behav- 
ior for regular and chaotic orbits. On the other hand, it is clear that following 
a great number of deviation vectors requires considerably more computation 
time. It turns out, however, that, if chaos occupies a 'large' portion of phase 
space, a well-tailored application of GALI^, with 2 < k < N, can significantly 
diminish the CPU time required for the detailed 'charting' of phase space, 
compared with SALI [k = 2), as we demonstrated on specific examples in 
section 15.21 (see figure [8]). 

Although the results presented in this paper were obtained for degree of 
freedom Hamiltonian systems, it is easy to see that they also apply to 2A^- 
dimensional symplectic maps. So, equations (j^3l) and flHOj) which describe the 
behavior of GALI^, with 2 < k < 2N, for chaotic and regular orbits respec- 
tively are expected to hold in that case also. One remark is in order, however: 
In the case of A^ = 1, i. e. for 2D maps, the first condition of equation flHUl) 
cannot be fulfilled. Thus, for regular orbits of 2D maps, any 2 initially inde- 
pendent deviation vectors will become aligned in the direction tangent to the 
corresponding invariant curve and GALI2 will tend to zero following a power 
law of the form GALI2 oc l/t"^. This behavior is already known in the literature 
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A Wedge product 



Following an introduction to the theory of wedge products as presented in 
textbooks such as [52], let us consider an M-dimensional vector space V over 
the field of real numbers M. The exterior algebra of V is denoted by A(V^) 
and its multiplication, known as the wedge product or the exterior product, 
is written as A. The wedge product is associative: 

{u Av) Aw = u A {v Aw) (A.l) 



for u^v^weV and bilinear 



(ciM + C2V) Aw = Ci{u Aw) + C2{v A w) , 

W A (CiM + C2V) = Ci{w Au) + C2{w A v) (A. 2) 

for u,v,w E V and Ci, C2 G M. The wedge product is also alternating on V 
uAu = (A.3) 

for all vectors u E V. Thus we have that 

uAv = -vAu (A.4) 

for all vectors u,v E V and 

■Ui A ^2 A ■ ■ ■ A Mfc = (A. 5) 



whenever ui,U2, . . . ,Uk £ V are linearly dependent. 

Elements of the form Ui A U2 A ■ ■ ■ A Uk with ui,U2, . . . ,Uk G V are called 
fc-vectors. The subspace of A(V^) generated by all /c-vectors is called the fc-th 
exterior power of V and denoted by A^{V). The exterior algebra A{V) can be 
written as the direct sum of each of the k-th powers of V: 

M 

A{V) = A\V) = A°(r) © A\V) ®A\V)®---® A^\V) (A.6) 

A;=0 



where A%V) = M and A\V) = V. 

Let {ei, 62, ... , Cm} be an orthonormal basis of V, i. e. Cj, i = 1, 2, . . . , M are 
linearly independent vectors of unit magnitude and 

ii ■ ij = 6ij (A. 7) 
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where (■) denotes the inner product in V and 
, 1 for i = 7 

S.J = { . (A.8) 

for i ^ j 



It can be easily seen that the set 

{ei, A A ■ ■ ■ A Ci, | 1 < ii < i2 < ■ ■ ■ < ik < M} (A.9) 

is a basis of A'^iV) since any wedge product of the form ui A U2 A ■ ■ ■ A Uk 
can be written as a hnear combination of the fc- vectors of equation (IA.9p . 



This is true because every vector Ui, i = 1,2, ... ,k can be written as a hnear 
combination of the basis vectors e,, i = 1,2, . . . , M and using the bihnearity 
of the wedge product this can be expanded to a hnear combination of wedge 
products of those basis vectors. Any wedge product in which the same basis 
vector appears more than once is zero, while any wedge product in which the 
basis vectors do not appear in the proper order can be reordered, changing 
the sign whenever two basis vectors change places. The dimension of A''{V) is 
equal to the binomial coefficient 



dimA'=(r) 



and thus the dimension of A{V) is equal to the sum of the binomial coefficients 



dimA(\/) = ^1^1=2*^. (A.ll) 



M 

k=0 \ k 



The coefficients of a /c-vector -Ui A'U2 A - ■ -Auk are the minors of the matrix that 
describes the vectors Ut, i = 1,2, k in terms of the basis e., i = 1,2, . . . , M. 
Let us write these relations in matrix form 



Ui 




Uii 


Ul2 ■ 


■ UiM 




ei 




ei 


U2 




U21 


U22 ■ 


■ U2M 






= C 


62 


Uk_ 




Ukl 


Uk2 ■ 


■ UkM 




eu _ 




_ 



(A.12) 



C being the matrix of the coefficients of vectors Ui,i = 1,2, ... ,k with respect 
to the orthonormal basis e,, i = 1,2, ...,M and Uij, i = 1,2, . . . , k, j = 
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1,2, ... ,M being real numbers. Then the wedge product -ui A ^2 A ■ ■ ■ A ma; is 
defined by 



Ui A U2 A ■ ■ ■ A Uk = 

l<il<i2<-<ife<M 



Uii^ Uii^ ■ ■ ■ Uii^ 



A A ■ ■ ■ A Ci^ (A 



where the sum is performed over all possible combinations of k indices out of 
the M total indices. So the coefficient of a particular fc- vector e^^ Ae^j A ■ ■ ■ ACj^. 
is the determinant of the k x k submatrix of the k x M matrix of coefficients 



appearing in equation (IA.12P formed by its ii, 12, ■ ■ ■, ik columns. 



B The relation between GALI2 and SALI 



Proposition 1 We consider a 2N -dimensional vector space over the field of 
real numbers M, which has the usual Euclidean norm and is spanned by the 
orthonormal basis {ii, 62, ... , e2N}- We also consider two unit vectors wi, W2 
in this space so that 



2N 2N 



(B.l) 



i=l 



i=l 



and 



2N 



2N 



w. 



2i 



(B.2) 



i=l 



i=l 



Let us now define the 2-vector Wi A W2 from equation liA.13\) and its norm 
from equation [T^) . Under these assumptions the following holds: 



\Wi A W2\ 



IWi - li)2 • + W2\ 



(B.3) 



Proof. Expanding the right hand side of equation (]B.3P we have: 



A 



\UJ1 + UJ2\\ ■ \\Wi - W2\ 



I Ei=l [Wli - U)2i) ■ Ei=l [Wli + W2i 
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1 

4 



\i=l 



2N 



27V 



'27V 



27V 



IN 



1=1 



( II + I] W^l - 2 ^ WiiW2i ) ■ (Y.^li+J2^l + ^Y1 ^li^2i 
i=l / \i=l i=l i=l / 

2N \ / 2N \ / 2Af \ ^ 

1 - ^ WliW2ij ■ y- + Yl ^li^2i j = 1 - [Yl ^li^2i j ^ 

/27V \ 
A=l- \J2 ^li'^li + 2 XI WliW2iWljW2j , (B.4) 



ii=l 



where we made use of f lB.2p . On the other hand, using equation f|T8l) we get 
for the left hand side of equation fIB.SP : 



B=\\w,AW2f = J2 



i<j 



W2i W2j 



Y{WiiW2j - WijW2i 



i<j 



^ = J2 ^li^lj + J2 ^lj^2i - 2 X WuW2iWijW2j. 

i<j i<j i<j 



(B.5) 



The first two sums of equation (IB.SP contain all the possible products of the 
coordinates of the two vectors except the ones corresponding to equal indices, 
i = j. So the quantity B can be written as follows: 



^ = Y1 ^li^lj - 2 X WiiW2iWijW2j 
ij^j i<j 

27V / 2N 

J2 w\-wlj + X W\(W\ - ^l^li + 2 I] WuW2iWxjW2j 

i^j i=l \i=l i<j 



(B.6) 



Now, the first two sums contain all the possible products between the coordi- 
nates of the two vectors and so B takes the form: 



2N 27V / 2N 

B = YY1 ^li^lj - Y ^li^li + 2 51 Wx^2iWxjW2j 
i=l j=l \i=l i<j 

27V 27V / 27V 

= Y^U-Y.^li- [Y. ^li^li + 2 I] Wx^2iWxjW2j 
\i=\ i<j 



i=l i=l 



2N 



B =1- \Y wliWli + 2 X WiiW2iWijW2 

i<j 



'3 ' 



.2=1 



(B.7) 



where we used again (]B.2|) . Comparing equations (1B.4P and (IB. 70 we see that 
both sides of equation (IB.SP are equal and so the proof of proposition [1] is 
complete. ■ 
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Using equation (1B.3P as well as the definitions of SALT ([6]) and GALI2 (1191) 
we conclude that the precise relation between the two indices is 

GALI. = SALI . + ^g^g^ 

So, the two indices are proportional to each other 

GALI2 oc SALI, (B.9) 

since the quantity m = max{||wi + W2II , \\wi — W2\\} lies in the interval m G 
[-\/2, 2]. In particular, in the case of chaotic orbits m — 2 as SALI 
and eventually GALI2 also vanishes, while in the case of regular motion m 
fluctuates around non-zero values in the above interval [a/2,2). 

From the above discussion we conclude that SALI is essentially equivalent to 
GALI2. In practice, however, since the computation of GALI2 according to 
equation (ITSi) ioi k = 2, requires the evaluation of several 2x2 determinants, 
it is more convenient to compute SALI in its place, by performing the simpler 
computation of equation IQ. 
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